## This code include scripts for Appendix Tables and Figures
# Table A18
# Figure A7


rm(list = ls())

library(readr)
library(ggplot2)
library(stargazer)
library(lme4)
library(dplyr)
library(haven)
library(texreg)
library(tidyr)

load("Final_JPR/data/paperdata.RData")
load("Final_JPR/data/paperdata_nokin.RData")
data <- paperdata %>%
       filter(tek_count>0)

source("Final_JPR/Code/coefplot.R")
source("Final_JPR/Code/marginal_effects.R")
source("Final_JPR/Code/setup.R")




data <- paperdata %>%
      plyr::mutate(tek_dum = ifelse(tek_count > 0, 1, 0))
table(data$tek_dum)
paperdata_nokin <- paperdata_nokin %>% dplyr::mutate(tek_dum = ifelse(tek_count > 0, 1, 0))
table(paperdata_nokin$tek_dum)

onset_logit <-  glm(onset_do_flag ~ lineq_total*tek_dum + lineq_total +
                           status_excl + lsize + family_warhistdummy + family_downgraded2 + 
                           ctygdppc_log + ctypop_log + peaceyears+ peaceyears_sq+ peaceyears_cub,
                    family = binomial(link = "logit"),
                    data = paperdata_nokin)


onset_fixed <-  lm(onset_do_flag ~ lineq_total*tek_dum + lineq_total  +
                          status_excl + lsize + family_warhistdummy + family_downgraded2 + 
                          ctygdppc_log + ctypop_log + factor(gwgroupid) + factor(year),
                   data = paperdata_nokin)

onset_mixed <-  glmer(onset_do_flag ~ lineq_total*tek_dum + lineq_total + lineq_total_sq + 
                     status_excl + lsize + family_warhistdummy + family_downgraded2 + 
                     ctygdppc_log + ctypop_log + (1 |year)+ (1 | countries_gwid),
              family = binomial(link = "logit"),
              data = paperdata_nokin)



incidence_logit <- glm(incidence_flag ~ lineq_total*tek_dum  + lineq_total + 
                              status_excl + lsize + family_warhistdummy + family_downgraded2 + 
                              ctygdppc_log + ctypop_log + peaceyears+ peaceyears_sq+ peaceyears_cub,
                       data = paperdata_nokin )


incidence_fixed <- lm(incidence_flag ~ lineq_total*tek_dum + lineq_total +
                             status_excl + lsize + family_warhistdummy + family_downgraded2 + 
                             ctygdppc_log + ctypop_log + factor(gwgroupid) + factor(year),
                      data = paperdata_nokin )


incidence_mixed <- glmer(incidence_flag ~ lineq_total*tek_dum + lineq_total +
                          status_excl + lsize + family_warhistdummy + family_downgraded2 + 
                          ctygdppc_log + ctypop_log + (1 |year)+ (1 | countries_gwid),
                   family = binomial(link = "logit"),
                   data = paperdata_nokin )




screenreg(list(onset_logit,onset_fixed,onset_mixed, incidence_logit,incidence_fixed,incidence_mixed),
          stars = c(0.01, 0.05, 0.1),
          caption = "Robustness Check: Heterogeneity analysis of Horizontal inequality Conditional on TEK (1992-2020)",
          caption.above = TRUE,
          label = "tab:tab11",
          scalebox = 0.7,
          use.packages = FALSE,
          custom.header = list("DV: Conflict Onset" = 1:3,
                               "DV: Conflict Prevalence" = 4:6),
          custom.model.names = c("M1:logit", "M2:fixed", "M3:mixed", 
                                 "M4:logit", "M5:fixed", "M6:mixed"),
          custom.coef.map = list("lineq_total" = "Horizontal inequality",
                                 "lineq_total:tek_dum" = "Horizontal inequality $\\times$ TEK(dummy)",
                                 "tek_dum" = "TEK(dummy)",
                                 # "lineq_total_sq" =  "Horizontal inequality $^2$",
                                 "status_excl" =  "Status excluded",  
                                 "lsize " =  "Relative group size",
                                 "family_warhistdummy" =  "Previous rebellions",
                                 "family_downgraded2" = "Status downgraded",
                                 "ctygdppc_log" = "Ln(Country GDP per capita)",
                                 "ctypop_log" = "Ln(Country population)",
                                 "peaceyears" = "Peace Year",
                                 "peaceyears_sq" =  "Peace Year $^2$",
                                 "peaceyears_cub" = "Peace Year $^2$",
                                 "(Intercept)" = "Intercept"),
          custom.gof.rows = list("Group Fixed Effects" = c("NO", "YES", "NO", "NO", "YES", "NO"),
                                 "Year Fixed Effects" = c("NO", "YES", "NO", "NO", "YES", "NO")),
          digits = 3)


# Table A18
texreg(list(onset_logit,onset_fixed,onset_mixed, incidence_logit,incidence_fixed,incidence_mixed), 
       file = "Final_JPR/tables/appendix_tables/Table_A18.tex",
       stars = c(0.01, 0.05, 0.1),
       caption = "Robustness Check: Heterogeneity analysis of Horizontal inequality Conditional on TEK (1992-2020)",
       caption.above = TRUE,
       label = "tab:tab11",
       scalebox = 0.7,
       use.packages = FALSE,
       custom.header = list("DV: Conflict Onset" = 1:3,
                            "DV: Conflict Incidence" = 4:6),
       custom.model.names = c("M1:logit", "M2:fixed", "M3:mixed", 
                              "M4:logit", "M5:fixed", "M6:mixed"),
       custom.coef.map = list("lineq_total" = "Horizontal inequality",
                              "lineq_total:tek_dum" = "Horizontal inequality $\\times$ TEK(dummy)",
                              "tek_dum" = "TEK(dummy)",
                             # "lineq_total_sq" =  "Horizontal inequality $^2$",
                              "status_excl" =  "Status excluded",  
                              "lsize " =  "Relative group size",
                              "family_warhistdummy" =  "Previous rebellions",
                              "family_downgraded2" = "Status downgraded",
                              "ctygdppc_log" = "Ln(Country GDP per capita)",
                              "ctypop_log" = "Ln(Country population)",
                              "peaceyears" = "Peace Year",
                              "peaceyears_sq" =  "Peace Year $^2$",
                              "peaceyears_cub" = "Peace Year $^2$",
                              "(Intercept)" = "Intercept"),
       custom.gof.rows = list("Group Fixed Effects" = c("NO", "YES", "NO", "NO", "YES", "NO"),
                              "Year Fixed Effects" = c("NO", "YES", "NO", "NO", "YES", "NO")),
       digits = 3)



#Figure A7


m1 <- plot_interEffect_binary(ModelResults =incidence_fixed, bvarname1 = "tek_dum",
                              cvarname2 = "lineq_total", clusterid = "gwgroupid", data = paperdata_nokin,
                              cval1 = -14, cval2 = 2.5, dv = "incidence_flag")

p = ggplot(m1) + 
       geom_hline(aes(yintercept=0), linetype=2, color = "black") + 
       geom_pointrange(aes(x = type, y =mean, ymin =lo,
                           ymax =hi, shape = IV),
                       lwd = 1, position = position_dodge(width = 1/2),
                       fill = "WHITE")+
       
       #geom_point(size=4,  position = position_dodge(width = 1/2)) + 
       geom_linerange(aes(x = type, ymin=lo, ymax=hi),
                      alpha = 1, size = 1,
                      position = position_dodge(width = 1/2)) + 
       theme_bw() + scale_x_discrete(labels = c("Without TEK","With TEK"))+
       xlab("") + ylab("Average Marginal Effects")+
       theme(legend.position="none",
             legend.title=element_blank(),
             axis.text = element_text(size=14),
             text = element_text(size=14),
             plot.title = element_text(hjust = .5, size = 12, face = "bold"))+
       guides(color = guide_legend(nrow = 2))+
       ggtitle("Marginal Effects of Horizontal Inequality Conditional on the Presence of TEK")


ggsave("Final_JPR/Figures/appendix_figures/fig_A7.pdf", width = 7, height = 4.2, units='in', dpi=600)
